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Expansion of a Bose-Einstein condensate (BEC) is studied, in the presence of 
a random potential. The expansion is controlled by a single parameter, (/j,T e ff/h), 
where [i is the chemical potential, prior to the release of the BEC from the trap, 
and T e jf is a transport relaxation time which characterizes the strength of the 
disorder. Repulsive interactions (nonlinearity) facilitate transport and can lead to 
diffusive spreading of the condensate which, in the absence of interactions, would 
have remained localized in the vicinity of its initial location. 

PACS numbers: 03.75Kk, 71.55 Jv 



There is a number of recent experiments on the interplay between nonlinearity and 
disorder [1, 2, 3, 4, 5, 6]. In [1, 2, 3, 4] the expansion of a BEC, in the presence of a 
random potential, has been studied. In [6] the authors study propagation of a light beam 
through a dielectric sample, with a random refraction index, under linear and nonlinear 
conditions. Under certain conditions both types of experiments can be described with the 
Gross-Pitaevskii equation. For the BEC problem the necessary conditions are specified, 
e.g., in [7]. For optical waves the Gross-Pitaevskii equation emerges within the paraxial 
approximation (see [6] and references therein). The equation is: 
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where m is the atom mass, g = Anah 2 jm is the interaction parameter related to the 
scattering length a, V(r) is the random potential and the wave function is normalized to 
the total number of atoms, 



We have written the equation in the BEC context and will use the appropriate language 
throughout the paper (a brief translation into the language of optics will be made at the 
end). The random potential can be generated optically (a laser speckle pattern, [8]) or it 
can be due to some foreign (impurity) atoms, with which the BEC interacts [9]. Eq. (1) 
is applicable in any spatial dimension d = 1,2,3. In the 2d case, to which most of this 
paper is devoted, the condensate can move only in the (x, |/)-plane, whereas the motion 
in the z-direction is frozen. For this geometry r represents the in-plane coordinate and 
A is the corresponding 2d Laplacian. The interaction parameter for the 2d problem is 
g = (87r) 1 / 2 (/i 2 a/ma z ), where a z is the oscillator length in the z-direction [7]. We assume 
positive g, i.e., repulsive interactions. Eq.(l) has to be solved with the initial condition 



where F(r) represents the initial shape of the condensate, i.e., the equilibrium shape in 
the trap prior to its release. It is assumed that at t — the trap is suddenly switched off 
and the expansion, described by Eq. (1), starts. 

We will be interested in strong nonlinearity in the sense that, initially, almost the 
entire energy of the condensate is due to interactions. In this case, and in the absence 
of the random potential (V(r) = 0), the expansion is known to occur in two stages 
[10, 11, 7]: a rapid "explosion", followed by expansion linear in time. The explosive 
stage takes time of order t = 1/uj, where u> is the frequency of the parabolic trap, from 
which the BEC is released. We assume that the trap is isotropic (for the 2c? case the 
isotropy pertains, of course, only to the (x,y)- plane but not to the z-direction). During 
the time to the condensate, roughly, doubles in size and most of its interaction energy 
gets converted into the flow energy. The second stage corresponds to times larger than to, 
when the non-linearity becomes weak and the expansion is accurately described by the 
linear Schrodinger equation. This observation suggests a simple, approximate approach to 
the problem of expansion in the presence of disorder: If the disorder is sufficiently weak 
(see the criterion below), then the first stage of the expansion is dominated by strong 
nonlinearity, while the disorder can be neglected. In the second stage the non-linearity 
can be neglected, so that one is left with the linear Schrodinger equation, in the presence 
of disorder. This idea has been proposed in [12], which appeared when the present work 
was in progress (in [12] the Id case was considered, whereas we are interested mostly in 
the 2d and 3d case). Thus, the evolution of the BEC at the second stage of the expansion 
amounts to a diffusion process in the weakly disordered potential. 

At the end of the first stage of the expansion, i.e., at time of order to, the condensate 
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can be roughly described by a wave packet [10, 11, 7] 

V(r,t ) ~ F(r)exp(ir 2 /al), (4) 

where a = (h/muo) 1 / 2 is the oscillator size of the (in-plane) trap. The function F(r), on 
a qualitative level, can be taken the same as the initial wave function in the trap: indeed, 
during the time to the condensate has only expanded by a factor of 2 or so, preserving 
its shape of the inverted parabola. The essential thing that happens during the time t 
is that the wave function acquires a phase, which accounts for the flow energy. Since 
the velocity of the condensate is proportional to the gradient of the phase in Eq.(4), it 
follows that the velocity is proportional to r, which is the main feature of the ballistically 
expanding condensate. The wave function in Eq.(4) serves as the initial condition for the 
second stage of expansion, controlled by the linear Schrodinger equation, in the random 
potential V(r). To facilitate the calculations below, we shall approximate the initial func- 
tion by a Gaussian, F(r) = Aexp(—r 2 /2Rl), where A = (N/ttRI) 1 ^ 2 is a normalization 
factor and Rq characterizes the initial size of the condensate. This replacement of the 
inverted parabola by a Gaussian does not affect the evolution of the wave packet (4) in 
any significant way. 

A clear separation into the two stages can be made only if the disorder is sufficiently 
weak [12]. This means that during the initial, "explosive" stage there is no scattering of the 
condensate on the random potential, i.e., this initial stage is ballistic. This requirement 
is satisfied if the transport mean free path / is larger than 2ir/k, where k is the wave 
number of the wave (particle) in question. The wave packet in Eq.(4) contains a broad 
range of fc's, with a Gaussian tail for k > k , where k = Rq/o^. The value of I separates 
the spatial Fourier components of the packet into two groups. The components with 
kl » 1 will evolve by a diffusion process, whereas those with kl < 1 will be localized by 
disorder and, thus, will stay in the vicinity of the initial location of the packet. (In 2d 
also the large-A; components will eventually get localized, but for those the localization 
length is exponentially large so that there is a lot of room for the diffusive propagation). 
In order to clearly observe the diffusive behavior one must take / 3> k^ 1 . This is a rather 
weak condition, because k^ 1 is of the order of the healing length £ which is the smallest 
macroscopic length in the problem. (Note that we assume that the bottom of the trap, 
from which the condensate was released, coincides with the average value of the random 
potential; if the bottom were raised sufficiently with respect to that value, then even 
the small- A; components would have enough energy to diffuse away.) To ensure that the 
great majority of the Fourier components would evolve by diffusion we choose, somewhat 
arbitrarily, / i? , so that the hierarchy of length scales is 

£ < a < R < l (5) 

Note that, since Rq ~ fi/muj 2 and £ 2 ~ h 2 /mfi, the ratios £/ao and a®/ Rq are of the same 
order, (hu/fi) 1 ^ 2 , where /j, is the chemical potential of the BEC before it was released from 
the trap. 
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Let us now present the calculation of the time evolution of a packet, whose initial 
shape is defined in (4), with F(r) being a Gaussian, i.e., 

tf(r,i = 0) = (N/nRl) 1/2 exp{-r 2 /2R 2 )exp(ir 2 /2al), (6) 

where to in the argument of \I/ has been replaced by 0. This is because we are studying now 
the time evolution starting from to (and on time scale much larger than to). Decomposing 
ty(r, t = 0) into Fourier components, ^f(k), we obtain 

(*0 = (N/irR 2 ) 1/2 b 2 exp(-b 2 k 2 /2), (7) 

where the complex length b is defined as 

1 



The quantity | SS/(k) \ 2 describes the density of the particles in k-space. More precisely, 
dn(k) =| ^f(k) | 2 d 2 k is the number of particles in the interval d 2 k, or in the interval d 2 Vk 
in the velocity space (vk = Tik/m). The distribution drik, emerging from the "explosion", 
will evolve by diffusion, due to scattering from the random potential V(r). One must keep 
in mind that the diffusion coefficient, D(k) = ^v(k)l(k) ( as well as the mean free path 
I), depend on the velocity, so that the group dn(k) of particles, after time t, will spread 
into a cloud 

This expression has to be integrated over k, to obtain the BEC shape, | \P(r, t) \ 2 , at 
time t. To perform this integration we have to specify the dependence of /, or of the 
mean transport time r, on k. For a Gaussian white noise potential, r is known to be 
inverse proportional to the density of states v{h) [13]. Since, in 2d, v = const, we have 
D{k) = \v 2 (k)r = |^-t, with r being a constant specifying the strength of the random 
potential. Integration over k of the expression (9) yields 

'^i 2 ^*^ (10) 

where K is the zeroth order modified Hankel function and D is the diffusion coeffi- 
cient corresponding to k — k , i.e., D = \v^t = |^-^-t. The mean square size of the 
condensate grows as 2D t. 

Eq.(10) describes the shape of the condensate in the diffusion regime. It is quite re- 
markable that this shape depends on the single quantity Do- It is instructive to rewrite 
Do in terms of the chemical potential fi of the BEC in the trap, prior to its expansion. 
Using kg = i?o/ a o = 2mfi/h 2 , we obtain D = fir/m. Thus, the shape (for a given m) 
is determined by a single dimensionless parameter /it/H. The situation turns out to be 
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analogous to a disordered Fermi system, where (Epr/h) ~ kpl is the only parameter 
which determines the transport properties of a disordered metal, with Fermi energy Ep 
[13]. The necessary condition for the diffusive regime described above is that the param- 
eter (fir/h) ^> 1. When the disorder increases and fir/h approaches unity the diffusion is 
inhibited and the BEC gets localized. (In 2d, even for (/rr/fr) 3> 1 spreading by diffusion 
will eventually stop and localization is expected to set in. The localization length, how- 
ever, is exponentially large, L[ oc ~ exp(fir/h).) Instead of changing the disorder, one can 
tune the interaction strength, i.e., the nonlinearity (for instance, using a Feshbach reso- 
nance), and observe a crossover from localization (weak interactions) to diffusion (strong 
interaction). 

In the derivation of Eq.(lO) we have ignored the fact that the small-/c components of 
the initial wave packet (6) cannot propagate by diffusion, as was already discussed above. 
To account for this fact one should integrate the expression (9) with a lower cutoff, k min . 
The cutoff is determined by the condition kl(k) ~ 1 , i.e., k min ~ (m/hr) 1 ^ 2 . Although 
it is not possible then to derive a close analytic expression for | \I/(r, t) | 2 , the qualitative 
picture is clear: most of the wave packet will evolve according to Eq.(lO), but a small 
portion will remain localized near the origin. Furthermore, we have neglected the weak 
localization effects, which influence the large-A; components and renormalize (by a small 
amount) the corresponding diffusion coefficient D{k). 

Our calculation can be extended to three dimensions. Eq.(9) remains the same, with 
d 3 k instead of d 2 k and with 

I ^j(k) \ 2 = - e~ k2/k o (11) 

1 [)l (W^) 3 ' 1 j 

The main difference with the two-dimensional case is that in 3d the density of states, v(k), 
is proportional to k. This leads to a linear dependence of the diffusion coefficient on k 
(instead of the square dependence, in 2d). The linear dependence complicates integration 
over k, and we were not able to obtain a close expression for | ^(r, t) \ 2 . The essential 
thing, however, is that the spread of | \&(r, i) \ 2 is still determined by the single quantity, 
Do = ^v l = ^VqTq, where r = r(/c ) (note that, unlike in 2d, r now depends on k). Thus, 
again, the behavior of the condensate is controlled by the single parameter (fj,r /h), and 
diffusion is possible only if this parameter is large. We will not dwell here on the diffusion 
regime but rather mention the interesting possibility of an interaction-induced Anderson 
transition. The point is that, in 3d, there is a critical value, c, such that for fiTo/h smaller 
than c the BEC remains localized. Therefore, by increasing the interaction strength, one 
can transfer the condensate from the localized regime into the extended one. 

Let us explain this crossover in more detail. For conceptual simplicity it is useful not 
to create disorder in the trap itself but only somewhat away. Then, in the absence of 
interactions /i would be small, of order hu, and the condensate released from the trap 
would not be able to penetrate the disordered region (we assume that Tiuj is smaller 
than the mobility edge E c ). With the increase of interactions /i, as well as the size of 
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the condensate Rq, increase and the wave function, after the first (ballistic) stage of 
expansion, develops oscillations (see Eq.(4)). The number of oscillations grows with /i, so 
that for larger \i the wave packet in Eq. (4) will contain more high energy components. The 
components with energy larger than E c will be able to penetrate the disordered region. 
Thus, under increase of interactions, larger and larger portions of the condensate will be 
transferred "over the mobility edge" and will diffuse away. 

The specific dependence of I (or r) on k is determined by the properties of the random 
potential. In our calculation a Gaussian, white noise potential has been taken. For a 
correlated potential the dependence of I on k will change (see, for instance, [14, 15]). Let 
us return to the 2d case and assume a Gaussian random potential, with zero mean, and 
with a correlation function < V(r)V(r r ) >= V^exipi— | r — r' | 2 /R^), where Vq and R c 
characterize, respectively, strength and range of the random potential. For this case a 
close expression for the transport mean free path / can be obtained [14]: 



where the function G(x) is given in terms of the modified Bessel functions as G(x) = 
exp(— x)[Io(x) — Ii(x)]. For kR c <C 1 (the white noise limit) I is proportional to k. In the 
opposite limit of a smooth random potential (kR c ^> 1) I is proportional to k A . D{k) in 
(9) becomes a complicated function of k and no closed expression for | \I>(r, t) | 2 can be 
obtained. The detailed shape of the expanding cloud is changed, as compared with the 
white noise case, Eq.(lO). The overall features, however, are preserved. For instance, the 
mean square size of the cloud is still proportional to D e fft, where the effective diffusion 
coefficient depends now on R c (in addition to V and k ). Thus, the rough picture of 
the expansion (but not the fine details) is universal and depends on the single quantity, 
D e ff. This statement pertains to short-range correlations, rapidly decaying beyond the 
distance R c (the correlation radius). It has been emphasized already in the original work 
of Anderson [16] that long range correlations in the randomness change the transport 
picture completely. The "one-dimensional" speckles, designed in [1, 2, 3, 4], do exhibit 
long range correlations, which has a significant effect on spreading of the Id condensate 
[12]. It is not clear how long range correlations will affect the diffusive spreading in 2 or 
3 dimensions. 

Let us briefly comment on the Id case. There is some numerical work on the subject 
[17, 18] claiming that even weak non-linearity (although still above a certain threshold) 
can lead to complete derealization. We should emphasize that weak nonlinearity effects 
are not addressed in the present work. We assumed strong nonlinearity, initially, and no 
nonlinearity after the first, explosive stage of the expansion, when almost all interaction 
energy had already been converted into flow. More recent numerics [19, 20] concentrated 
on strong nonlinearity, as appropriate for the experimental conditions [1, 2, 3, 4], and 
identified various regimes of Id localization, depending on the degree of disorder. In Id 
the localization length, L; oc , is of the same order as the mean free path I, so that there is no 




(12) 
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room for diffusion and our results do not apply. However, the classification of the disorder 
as "weak" or "strong", depending on the value of the parameter UqI ~ (fj,To/h), remains 
meaningful. If this parameter is large, then the details of the randomness, such as the cor- 
relation radius of the random potential, are not important (single parameter scaling [21]). 
When this parameter approaches unity, one enters the regime of strong localization where 
details become important and single parameter scaling is violated [22]. For a correlated 
Gaussian potential, < V(x)V(x') >= lg 2 /(| x — x' \ /R c ), the inverse localization length, 
for a given wave number k, is of the order of [23] Lj^. ~ (mV /h 2 k) 2 f(kR c ), where f(kR c ) 
is the Fourier transform of some (short range) correlation function /(| x — x' | /R c ). 
Note that f(kR c ) has units of length, i.e., it can be written as R c g(kR c ). For kR c <C 1, 
g(kR c ) ~ 1 (the white noise limit). In this limit L\ oc is proportional to k 2 /VqR c . For 
kR c ^> 1, Li oc rapidly increases, due to the decrease of the function f(kR c ), so that the 
"weak disorder condition", kLi oc 3> 1, becomes less restrictive. 

As has been mentioned at the beginning, Eq.(l) describes also propagation in nonlinear 
optics, in the paraxial (scalar wave) approximation. A monochromatic beam propagates 
in the z-direction and spreads in the transverse, (x, y)-plane. To translate Eq.(l) into 
the language of optics, one should make the following replacements: t =>- z, h =>• 
(c/u),m =>- l,V(r) =>- (—Sn(x,y)/no), where ui is the frequency of the wave, c is its 
speed (in the medium), no is the average refraction index, and 5n is its fluctuating part 
(it is crucial that 5n does not depend on z). Fluctuations of the refraction index in the 
transverse plane can lead to "transverse localization" [24] of the wave. The authors of [6] 
have recently observed this phenomenon, as well as the interplay between nonlinearity and 
disorder. The self-focusing (self-defocusing) in optics corresponds to attractive (repulsive) 
interactions in the BEC. Thus, the results presented in this paper applies to optics as well, 
for a self-defocusing nonlinearity. 

In conclusion, the diffusive spreading of a BEC, in the presence of a random potential, 
has been studied. The spreading is controlled by a single parameter, which can be conve- 
niently written as (/iT e ff/h). Here r e // is some effective relaxation time which depends on 
the properties of the random potential, as well as on u (i.e., on k = ^2mji/Ti). Diffusive 
spreading is possible only if this parameter is large. There is a clear parallel with the 
single parameter scaling in disordered electronic systems. There are also essential differ- 
ences between the two problems. In disordered conductors (at low temperatures) only the 
electrons in the vicinity of the Fermi energy participate in transport, whereas spreading 
of the BEC involves a broad range of wave numbers k (we considered only the case when 
the center of mass of the condensate is fixed, i.e., the condensate spreads but does not 
move as a whole). Strong repulsive non-linearity, considered in this paper, has a large 
effect on the BEC spreading and serves as a strong "delocalizing factor". 

Useful discussions with E. Akkermans, A. Minguzzi, M. Segev and J. Steinhauer are 
gratefully acknowledged. The research was supported by a grant from the Israel Science 
Foundation. 
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